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We study the dynamical evolution of disk and halo globular clusters in the Milky 
Way using a series of Fokker-Planck calculations combined with parametric statistical 
models. Our sample of 113 clusters with velocity data is predicted to descend from 
an initial population of 250 clusters, implying more than a factor of two decrease in 
population size due to evolution. 

Approximately 200 of these clusters are in a halo component and 50 in a disk 
component. The estimated initial halo population follows a coreless i? -3 38 density 
profile in good agreement with current estimates for the distribution of halo field stars. 
The observed core in the present-day distribution of halo clusters results from the rapid 
evaporation of clusters in the inner regions of the Galaxy. The initial halo population 
is also predicted to have a radially biased orbit distribution in rough agreement with 
the observed kinematics of halo field stars. The isotropy of the present-day halo cluster 
distribution results from the evaporation of clusters on elongated orbits. Similarly, the 
initial disk component has a nearly isotropic initial distribution that becomes more 
tangentially biased with time. However, the inferred initial characteristics of the disk 
component do not match the kinematics of the rapidly rotating thin or thick disk 
stellar populations. These characteristics may be more indicative of the flattened halo 
component discussed by Zinn (1993). 

Detailed examination of cluster evolution confirms the importance of disk heating. 
Clusters on low-inclination orbits experience the strongest disk heating because of 
optimal matches in resonant frequencies. Disk heating on high- inclination orbits is 
weaker but still dominates over spheroidal heating. Evaporation times depend weakly 
on initial concentration, density and height of oscillation above the disk. 

Key words: globular clusters: general - galaxies: individual (Milky Way) - galaxies: 
star clusters 



1 INTRODUCTION 

Age estimates for globular clusters indicate that they formed 
early in the history of the Milky Way and represent 'fos- 
sil relics' of the proto-Galaxy (Larson 1990). Attempting 
to uncover this history, researchers have carefully exam- 
ined a range of properties of the present-day cluster sys- 
tem, paying particular attention to the cluster kinematic 
distribution (e.g. Zinn 1993), mass distribution (e.g. Harris 
& Pudritz 1994), metallicity distribution (e.g. Zinn 1985) 
and age distribution (e.g. Chaboyer, Demarque & Saraje- 
dini 1996). These investigations have provided evidence for 
both accreted and native components in the cluster system 
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(Searle & Zinn 1978) and correlations between kinematics 
and metallicity which may trace the collapse of the Galaxy 
(Zinn 1985; Armandroff 1989; Zinn 1993). 

Comparisons with other stellar populations have re- 
vealed subtleties in the process of Galaxy formation and 
evolution. For example, in the inner Galaxy, the globular 
cluster distribution is flatter than the distribution of halo 
field stars although their profiles match well at larger radii. 
In addition, the field star velocity ellipsoid has a strong a ra- 
dial bias in comparison to the approximately isotropic clus- 
ter velocity ellipsoid (e.g. Ostriker, Binney & Saha 1989). 
Conversely, Zinn (1985) and Armandroff (1989) have pre- 
sented convincing evidence for a high-metallicity disk cluster 
population which has broad similarities in kinematics, spa- 
tial distribution and metallicity with the stellar thick disk. 
Understanding the origin of these relationships will improve 
our picture of the primordial Milky Way. 
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At the same time, theoretical interest in globular clus- 
ter evolution has been motivated by the discovery that two- 
body relaxation would drive evolution on a time scale that 
is much less than the age of a typical cluster (Ambartsum- 
ian 1938; Spitzer 1940; Chandrasekhar 1942). Subsequent 
research provided an understanding of the gravothermal in- 
stability (Lynden-Bell & Wood 1968) and the phenomenon 
of core collapse (e.g. Cohn 1980). One of the basic conclu- 
sions of this work on cluster evolution is that relaxation 
inevitably leads to evaporation (e.g. Spitzer 1987). Addi- 
tional refinements to the picture of relaxation-driven evolu- 
tion have been required to account for a source of energy 
which halts core collapse (e.g. Henon 1961; Lee & Ostriker 
1987) and to include tidal influences which arise on a clus- 
ter's orbit in a parent galaxy (Ostriker, Spitzer & Chevalier 
1972; Chernoff, Kochanek & Shapiro 1986; Weinberg 1994; 
Gnedin & Ostriker 1996; Murali & Weinberg 1996, hereafter 
Paper I). 

Recent work on tidal influences has shown that evap- 
oration is accelerated by the interaction of a cluster with 
the tidal field produced by the halo and disk of the Galaxy 
(Gnedin & Ostriker 1996; Paper I). These studies show that 
depletion depends strongly on cluster mass and orbit in the 
Galaxy. This suggests that understanding the initial char- 
acteristics of clusters and their relationship to other stellar 
populations in the Galaxy requires a comprehensive descrip- 
tion of evolution since the time of formation. In a related 
study which supports this view, Murali & Weinberg (1996; 
hereafter Paper II), have demonstrated the importance of 
evolution in shaping the M87 globular cluster population 
and as a partial cause for the specific frequency conundrum 
in fundamental plane ellipticals. 

Motivated by these results in the present work, we 
investigate the degree to which evolution has shaped the 
Milky Way cluster population. Our approach is to predict 
the initial spatial and kinematic distributions of the glob- 
ular cluster system using the Fokker-Planck description of 
cluster evolution discussed in Paper I in combination with 
the parametric statistical framework employed in Paper II. 
The predictions describe the initial population of clusters 
which evolve quasi-statically through relaxation and tidal 
heating and indicate changes which dynamical evolution has 
wrought on the system as a whole. The results also provide a 
basis for understanding the primordial relationship of glob- 
ular clusters to other stellar populations. 

We first study the evolutionary behavior of clusters 
which inhabit the disk and halo of the Galaxy. The calcula- 
tions demonstrate the importance of disk heating on cluster 
evolution and quantify dependences on important internal 
and external parameters, including orbit in the Galaxy and 
cluster concentration. We also examine the behavior of in- 
ternal density profiles and mass spectra in evolving clusters. 

Having considered the detailed physical behavior, we 
examine properties of the full cluster population. We first 
characterize properties of the current cluster population and 
then predict its initial conditions using the data set com- 
piled by Gnedin & Ostriker (1996) and the three-space ve- 
locities derived by Cudworth (1993). The inferences are de- 
rived from both spherical and two-component disk+sphere 
models of the cluster distribution. While several analyses 
have shown the cluster system to be approximately spher- 
ically distributed (Chernoff & Djorgovski 1989; Thomas 



1989), other investigations show two components: a flat- 
tened, rapidly rotating high-metallicity component associ- 
ated with the Galactic disk and a spherically distributed 
low-metallicity component associated with the Galactic halo 
(Zinn 1985; Armandroff 1989). Further subdivisions may 
also exist (Zinn 1993; Zinn 1996). The choice of models re- 
flects the gross characteristics of the cluster system and al- 
lows us to compare the candidate distributions. The results 
predict significant differences in the initial and present-day 
cluster populations, indicating the role of evolution in shap- 
ing the present-day cluster system. Moreover, neither model 
is completely successful in describing the cluster population 
probably due to the combined effects of evolution and ob- 
scuration. 

The plan of the paper is as follows. In we summarize 
the approach and scenario used throughout the investiga- 
tion. The results are presented in §^|and include description 
of the physical behavior as a function of orbit, examination 
of the internal properties of evolving clusters and analysis 
and prediction of the initial conditions of the observed pop- 
ulation. Finally, §^ discusses the implications of the results. 
The appendices provide derivations of the models and a dis- 
cussion of the statistical procedure. 



2 SCENARIO AND INVESTIGATION 
2.1 Model populations 

Our fiducial population consists of clusters which formed in 
a single episode approximately HGyr ago, the lower limit 
on cluster ages estimated from current models of stellar evo- 
lution (Chaboyer 1995) . For older ages, the evolution in the 
properties of the cluster system is greater than the estimates 
derived below. Initial clusters are assigned Wo = 5 King 
mo del pr ofiles. Investigation of concentration dependence 
in §3.1.5 shows that evaporation times vary little with Wo. 
We assume that each cluster has a Salpeter IMF m - ^ with 
— 2.35 and lower mass limit mi —0.1 Mq. For this choice, 



stellar evolution would dominate for the first Gyr, roughly 
corresponding to the main sequence lifetime of a 2Mq A- 
star, which we choose as the upper mass limit, m u . Follow- 
ing the phase of strong stellar evolution, relaxation, external 
heating and, ultimately, core collapse heating would begin 
to drive cluster evolution. We define our zero-population at 
this epoch, approximately 10 Gyr in the past. 

In light of observational evidence, we adopt a two- 
component model of the cluster population consisting of 
flattened and spherical distributions. We employ the com- 
monly used terminology of 'disk' and 'halo' cluster to refer 
to members of these sub-populations. A 'disk' cluster is most 
often 'metal-rich' with disk kinematics while a 'halo' cluster 
is most often 'metal-poor' with halo kinematics but which 
may have high or low orbital inclination relative to the disk. 
A 'classic' halo orbit, however, is one of high inclination from 
the disk. 



2.2 Cluster evolution 

Following formation and early stellar evolution, cluster evo- 
lution is driven by relaxation, the tidal field and binary heat- 
ing of the core. As described in Papers I & II, the compe- 
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tition between relaxation and the tidal field of a galaxian 
spheroid is particularly important in determining a cluster's 
evolutionary time scale and survival history. The dominant 
effects of the spheroid on tidally-limited clusters were found 
to be heating on low-eccentricity orbits and tidal truncation 
on high-eccentricity orbits. 

In the Milky Way, the disk also contributes significantly 
to cluster evolution. For high inclination orbits, a number of 
investigations have shown that the compressional shock im- 
parted to a cluster during its passage through the disk will 
generally enhance the evaporation rate (Spitzer & Chevalier 
1972; Chernoff, Kochanek & Shapiro 1986; Weinberg 1994; 
Gnedin & Ostriker 1996). For orbits confined to the disk, 
oscillations of the cluster about the midplane transfer en- 
ergy through resonant stellar orbits. Below we examine the 
effect of the disk on clu ster evolution on both low- and high- 
inclination orbits (§3.1). Appendix ^ outlines the derivation 
of the heating rate for disk oscillations. 



2.3 Orbits 

For eccentric, high-inclination orbits, the orbital phase of 
disk passage varies due to the precession of the argument of 
perihelion with respect to the plane of the Galaxy because 
orbits are not closed in the logarithmic potential. To remove 
dependence on orbital phase, we assume that disk shocking 
occurs twice per orbital period at the average orbital radius 
of the cluster. For eccentric low-inclination orbits, we follow 
the same procedure and also assume that the vertical motion 
is separable from the radial and tangential motion in all 
orbits. This allows us to define approximate three-integral 
distribution func tion s (DFs) in terms of algebraic constants 
of the motion ( §|2.5| ). With these DFs, we can follow the 
evolution of the phase space distribution of gl obul ar clusters 
using a series of Fokker-Planck calculations (§2.8) 

Orbits in the spheroid are defined using the quantity 
k = J I J max (E), the angular momentum relative to maxi- 
mum for an orbit of energy E. The inclination angle, i, is 
defined with respect to the disk so that i = 0° defines an 
orbit in the plane of the disk. Oscillations through the disk 
are defined by their oscillation height in multiples of the disk 
scale height, zo- 



2.4 Tidal limitation 

Clusters on orbits highly inclined from the disk are tidally 
limited by the Galactic spheroid. While initial cluster den- 
sities may differ from the mean density required by peri- 
galactic tidal limitation, subsequent evolution during the 
first gigayear leads rapidly to tidal truncation or disruption. 
Clusters on orbits confined to the disk may posess limit- 
ing radii, R c , smaller than that implied by tidal limitation 
in the spheroid (e.g. note the discrepancy between the ob- 
served and predicted tidal radius of M71: Drukier, Fahlman 
& Richer 1993). This implies that their density is higher 
than that requi red fo r tidal limitation at given concentra- 
tion. However, 



3.1.3 



shows that evaporation times vary lit- 
tle with density for fixed mass. Therefore, in studying popu- 
lation evolution, we assume that low-inclination clusters are 
also tidally limited. The limiting radius R c is determined by 
the cluster mass, orbit and the ratio of cluster mean density 



to the mean density required for tidal limitation (see § 2. 1C 
and Paper I for further details). When clusters are tidally 
limited, we refer to the limiting radius as the tidal radius, 
Rt- 

2.5 Distribution functions 

We use parametric models to define distributions of cluster 
orbits and masses in both disk and spheroid populations. 
Mass distributions, u(M), are defined using either a power 
law (Harris & Pudritz 1994), a two-component power law or 
a Gaussian magnitude distribution (McLaughlin, Harris & 
Hanes 1994). The two-component power-law mass spectrum 
is continuous at the break at mass M cu t- 

Distributions of cluster orbits are defined using well- 
known models (Table |l|) which are generalized to provide 
kinematic and spati al d istributions for populations in the 
Galactic potential (§ |2.9[ ). Disk component DFs are defined 
as functions of Ed, the orbital energy associated with mo- 
tion in the disk plane, jf, the square of the z-component 
of the angular momentum and E z , the orbital energy as- 
sociated with motion perpendicular to the disk plane. The 
functional form is obtained using the Mestel disk with no 
net rotation (Binney & Tremaine 1987) combined with an 
isothermal vertical distribution. This provides a family of 
constant anisotropy, power-law surface density profiles with 
power- law index — (??d — 2q<i) of infinite range. The quantities 
r/d = Vc/ a d an d Id respectively define the velocity disper- 
sion, ad, in terms of the (constant) circular rotation velocity, 
v c , and degree of anisotropy of the distribution in the disk 
plane. The parameter rj z = $o/o"z defines the vertical scale 
height of the distribution in terms of the vertical velocity 
dispersion, a z , and the central potential of the disk, $o- Ap- 
pendix |c| discusses the model family in further detail. 

Spheroidal component distribution functions are de- 
fined as functions of E, the total energy of the orbital mo- 
tion, and J 2 , the square of the total angular momentum. 
We use two models. The first is obtained by adapting the 
Mestel DF to the spherical case which we refer to as the 
Mestel sphere. The second is the Eddington model (Aguilar, 
Hut & Ostriker 1988). The Mestel sphere provides a family 
of constant anisotropy power-law space density profiles with 
power- law index — (77 — 2q) of infinite range. The quantities 
77 = v'i/a' 2 and q respectively define the velocity dispersion, 
a, in terms of the (constant) circular rotation velocity, v c , 
and degree of anisotropy of the distribution as in the disk 
case. The Eddington sphere produces a family of variable 
anisotropy, power-law space density profiles with core radius 
R a and power-law index — 77 at large radii. Radial anisotropy 
becomes significant beyond the core radius. Appendices [Dj 
and [i] provide derivations and further discussion of these 
models. 

Complete orbit and mass distributions are given by joint 
distributions v(M) x /(/) where / denotes constants which 
define a disk or halo orbit. We assume that disk and halo 
components can have different mass spectra. In the two- 
component model, F denotes the fraction of the populat ion 
belonging to the spherical component. The results of §3.1 
show that disk cluster evolution varies little with oscillation 
height. Since the cluster population is not large, we therefore 
determine the parameter r\ z for the observed distribution, 
and keep it fixed when estimating the initial conditions in 
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order to minimize the number of parameters. For the disk 
distribution, we impose cutoffs in radius at Rd = 15 kpc 
(Wainscoat et al. 1992) and height at Z = 7.5 kpc, where Rd 
and Z are the radius in and height above the disk, respec- 
tively (see Appendix H for listing of coordinate notations). 



2.6 Statistical procedure 

A maximum likelihood method is used to estimate the pa- 
rameters in Table ^ from the observed cluster data. Our 
models can include 7 parameters for each cluster (mass, 3 
spatial coordinates and 3 velocity coordinates). However, 
not all quantities are available in most cases. To incorporate 
all available data, we use a likelihood technique for incom- 
plete data sets (Little & Rubin 1987; Stuart & Ord 1991). 

The likelihood is constructed in the usual manner as the 
joint probability of all observations given the model. How- 
ever, for clusters with unavailable phase space quantities, 
the probability is obtained by integrating the distribution 
function over all dimensions of unknown information. This 
defines the marginal probability of observing a particular 
cluster given the model. Typically, only the heliocentric ra- 
dial velocity is known, so we integrate over the tangential 
velocities relative to our vantage point to derive the marginal 
probability of observing a cluster at a given position with a 
given heliocentric radial velocity. Appendix |f] discusses the 
technique in more detail. 

Two types of fit are used to characteri ze th e data. Fits 
without the evolutionary calculations (c.f. §2.8) are used to 
derive the observed or present-day properties of the clus- 
ter sample. Fits based on the evolutionary calculations are 
used to derive the most likely initial conditions which pro- 
duce today's distribution. We first consider the present-day 
characteristics of the cluster system to serve as a guide to in- 
terpreting the evolutionary case. The listed uncertainties are 
1 — a variances computed under the assumption of normally 
distributed errors. 



as a shock on high-inclination orbits and as an average heat- 
ing rate on low-inclination, oscillatory orbits. The numerical 
procedure used on low-inclination orbits is the same used 
for orbits in the spheroid. On high-inclination orbits, we use 
an analogous, flux-conserving scheme to compute the total 
change in the DF due to the passage through the disk. 

Properties of the evolved cluster population are esti- 
mated from a grid of Fokker-Planck calculations performed 
over a range of cluster orbits. For the disk clusters, we use a 
4x4x5 grid in apogalactic radius, mass and k to sample the 
phase space. Apogalactica are taken in the range 2 < R a < 
8 kpc, masses in the range 10 s M Q < M < 5 x 10 6 M Q and 
orbits in the range 0.3 < k < 1.0. For the halo clusters, we 
use a 5 x 4 x 5 grid with 5 < R a < 15 kpc and the same range 
of mass and n. All clusters with R a = 0.8 kpc are assumed 
depleted and all clusters with energies equal to a circular 
orbit with R a = 20 kpc are assumed to be unevolved. For 
R a = 0.8 kpc, clusters on circular orbits either evaporate or 
decay into the nucleus by dynamical friction (Aguilar, Hut 
& Ostriker 1988) and clusters on eccentric orbits evaporate. 
For R a = 20 kpc, the evaporation timescale on a circular 
orbit is 70 Gyr for 10 5 Mq and is roughly the same for other 
orbits of equal energy. In most cases observed clusters with 
M < 10 s M Q initially had masses above 10 s Mq . In a few 
cases, extrapolation beyond the initially defined grid is re- 
quired. 

In § |3.l| , we find that low inclination halo clusters evolve 
more rapidly than high inclination clusters, although the 
differences are not extreme. However, in order to reduce 
computational expense, we neglect low-inclination orbits in 
the spherical component and assume that all orbits have 
high inclination. This assumption circumvents the roughly 
factor-of-four increase in halo phase space grid size required 
to sample the cylindrical geometry. As a result, the initially 
spherical distribution remains spherical and we underesti- 
mate the amount of evolution about the midplane of the 
disk. The consequences of this assumption are elaborated 
below. 



2.7 Data 

We use the data compiled in Gnedin & Ostriker (1996) which 
consists of 119 objects. Comparison with the Harris (1996) 
compilation shows no obvious systematic differences. We 
examine the distribution of clusters within 65 kpc having 
masses in the range 2.0 x 10 4 M Q < M < 2.75 x 10 6 M Q 
for a mass-to-light ratio of 3. This removes 6 clusters from 
the original sample, leaving a total of 113 clusters. Using 
the likelihood procedure outlined above, we also include the 
three-space velocities which now exist for about 20 clusters 
(Cudworth 1993). This is referred to as the augmented data 
set. 

2.8 Calculations 

We follow the evolution of individual clusters using the 
one-dimensional Fokker-Planck approximation (Cohn 1979) . 
The calculations include relaxation, external heating due to 
the time-varying tidal field of the disk and spheroid, and 
a phenomenological binary heating term (Lee et al. 1991). 
Spheroidal tidal heating is included using the implementa- 
tion discussed in Paper I. Disk tidal heating is implemented 



2.9 Galactic model 

We represent the spherical component of the Galaxy as a sin- 
gular isothermal sphere with Vb = 220 km s" 1 and the disk 
component as an exponential disk with radial scale length 
Ro = 3.5 kpc normalized to the disk central density in the 
solar neighborhood, p = 0.15M Q /pc 3 (Bahcall 1984). The 
vertical profile is taken to be Gaussian with a scale height 
of zq — 320 pc. The Gaussian was adopted instead of an 
exponential initially due to concern about the analyticity of 
the perturbation and its effect on adiabatic invariance. How- 
ever, comparisons between the two profiles show no strong 
differences in heating rate (Weinberg 1994). The solar ra- 
dius is taken to be Ro = 8.5 kpc. All measured radial veloc- 
ities are converted to velocities in the Galactic rest frame 
using an LSR velocity of 220 km s" 1 and a solar motion of 
(n, O, Z) = (-9, 12, 7) kms- 1 (Mihalas & Binney 1981). 

2.10 Parameterization of disk strength 

In order to parameterize the strength of the spheroid relative 
to the cluster in Paper I, we introduced the quantity M(x p ) 
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Table 1 . Functional forms of models 



Name 



Mass Models 
u(M) oc 



Parameters 



power law 

two-component 

power law 

Gaussian 



AI~ ai (M < M cut ) 
A/-«2(M > M C ut) 



012, Mcut 



-iV-V f/2«l. d y /dM 



cry 



Distribution Functions 
Name f(E, J 2 , E z ) <x Parameters 

Mestel disk e - E d/°l j^id e -E z /*l < % < oo 

< »7 Z < 00 

Mestel sphere e _B / ' 2 J 2 9 < 7) < 00 

— 1 < q < 00 

Eddington sphere e~ E l' y ' 1 e~ j2 / 2r l a ' 2 < q < 00 

< r a < 00 



which denotes the fraction of the cluster mass contained 
within the pericentric inner Lagrange point. Here, because 
of the the disk, it is convenient to define the following two 
ratios: \ = x p /R c denotes the ratio of the pericentric inner 
Lagrange point to the limiting radius of the cluster; and 
po,d{R) / pc denotes the ratio of the disk central density at 
a given radius R to the cluster mean density. This latter 
parameter defines the tidal amplitude of the disk relative to 
the cluster in the same way that M(x p ) is used to define 
the relative tidal amplitude of the spheroid. Using x> we 
can relate the cluster mean density, p c , to the mean density 
required for tidal limitation in the spheroid, p t , by p c /pt = 

x 3 . 



3 RESULTS 

3.1 Physical Behavior 

3.1.1 The importance of disk heating 

In Paper I, the evolution of a cluster of fixed mass, M(x p ) 
and k could be scaled to an orbit of any energy due to the 
scale-free nature of the tidal field. This allowed us to com- 
pare the orbital dependence of cluster evolution in several 
ways using the same calculations (c.f. Paper I, §3.1). Adding 
the disk destroys this scaling freedom because the quantity 
Po,d(Rd)/pc varies with the radius of disk crossing R d . In or- 
der to compare cluster evolution on different orbits, we show 
the mass remaining after 10 Gyr in tidally limited clusters 
on orbits of equal apocenter over a range of mass and eccen- 
tricity both with and without the disk (Tables |^ and ^|). 

For equal apocenter and fixed mass, the densities of 
tidally limited clusters increase with orbital eccentricity due 
to the decrease in Rt with increasing perigalactic angular 
frequency. Clusters on eccentric orbits therefore undergo the 



most rapid relaxation and have correspondingly large evapo- 
ration rates. As a result, for evolution in the spheroid alone, 
the remaining masses of 10 5 Mq clusters decrease monotoni- 
cally with eccentricity. Tidal heating also enhances mass loss 
rates on low-eccentricity orbits. The effect is noticeable in 
the remaining masses of more slowly relaxing 10 6 Mq clus- 
ters. At high eccentricity, relaxation still predominates, but 
at low eccentricity, heating becomes important. As a result 
there is a peak in remaining mass at intermediate eccentric- 
ity. 

The large tidal amplitude of the disk relative to the 
spheroid greatly enhances heating on low- and intermedi- 
ate eccentricity orbits in the disk+sphere calculations. High- 
mass low-eccentricity disk clusters lose equilibrium and dis- 
rupt. The low-mass counterparts do not disrupt but are 
rapidly driven to evaporation by strong tidal stripping. 
In addition, since the strength of the disk relative to the 
spheroid varies with radius, the importance of disk heating 
varies as a function of the radius of disk crossing. The rela- 
tive strength of disk heating is highest at about 8 kpc. As a 
result, the 1O 6 M0,k = 0.9 disk cluster at 8 kpc loses more 
mass than does its counterpart at 4 kpc (Table ^) . The tidal 
effect of the disk diminishes with increasing eccentricity be- 
cause cluster densities increase relative to the disk density 
at the crossing point. Disk heating is negligible at k = 0.3. 

Halo clusters exhibit the same overall tendencies as disk 
clusters, but heating rates are lower because resonances are 
concentrated at higher frequencies than resonances due to 
disk oscillations. For example, at 90° inclination, the lzo 
passage time scale is on the order of 1 Myr while the corre- 
sponding period of an oscillation with height greater than 
lzo is greater than 10 Myr for a disk cluster with Rd > 2 kpc. 
The location of the resonant frequency match is significant 
because the increase in binding energy at higher frequency 
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Table 2. Disk clusters: fraction of remaining mass after lOGyr 

10 s Mq 



Table 3. 



K 


1.0 


0.9 


0.6 


0.3 




11 


0.0 


0.0 


0.0 


0.0 


disk+sphere 




0.42 


0.35 


0.0 


0.0 


sphere 


8kpc 


0.0 


0.0 


0.05 


0.0 


disk+sphere 




0.77 


0.76 


0.40 


0.0 


sphere 


10 6 M Q 


K 


1.0 


0.9 


0.6 


0.3 




I kp 


0.0 


0.54 


0.74 


0.60 


disk+sphere 




0.77 


0.84 


0.89 


0.60 


sphere 


8kpc 


0.0 


0.25 


0.75 


0.86 


disk+sphere 




0.85 


0.91 


0.95 


0.86 


sphere 


for 52:0 oscillation height and Wo = 5 




ning mass 


after lOGyr 














10 5 M 






K = 


1.0 


0.9 


0.6 


0.3 




5 kpc 


0.14 


0.10 


0.0 


0.0 


disk+sphere 




0.58 


0.55 


0.0 


0.0 


sphere 


lOkpc 


0.74 


0.73 


0.42 


0.0 


disk+sphere 




0.82 


0.82 


0.59 


0.0 


sphere 


10 6 M Q 


K = 


1.0 


0.9 


0.6 


0.3 




5 kpc 


0.49 


0.56 


0.82 


0.72 


disk+sphere 




0.80 


0.86 


0.91 


0.73 


sphere 


10 kpc 


0.80 


0.82 


0.90 


0.90 


disk+sphere 




0.87 


0.93 


0.96 


0.90 


sphere 


for i = 9C 


and W = 


5 







reduces the response of a stellar orbit to a perturbation of 
fixed amplitude. 



3.1.2 Inclination dependence 

The importance of resonance location and frequency match- 
ing is also evident when comparing halo clusters on orbits of 
different inclination with respect to the disk. Table ^ com- 
pares cluster evolution on circular orbits and shows that 
mass loss increases as orbits become less inclined because 
resonances appear at lower frequency. At 5 kpc, the best 
match occurs for i — 30°, while at 10 kpc, the best match 
occurs at i — 15°. Radial differences arise because the disk 
passage time is fixed at lzo/v c sini while the cluster dynam- 

-—1/2 

ical time scale varies as p c as determined by the mean 
density of the spheroid through tidal limitation. The drop in 
cluster mean density at 10 kpc leads to more efficient heating 
at lower passage speed. 

The tendency for clusters confined to the disk and at 
low inclination to evolve more rapidly implies that an ini- 



tially spherical halo distribution will develop a vertically- 
dependent density profile with minimum density about the 
midplane of t he dis k. This effect is most pronounced at low 
eccentricity (j j3.1.1 ) because the high densities of tidally 
truncated clusters on highly eccentric orbits strongly re- 
duce the resonance amplitudes between stellar orbits and 
the time-varying disk tidal field. 



3.1.3 Density dependence 

The rapid disruption of tidally-limited low-eccentricity disk 
clusters implies that only those with high initial densities 
can survive for long periods after formation. Figure [l] shows 
survival, disruption and evaporation patterns in clusters on 
circular orbits with a range of initial densities, parameter- 
ized by the ratio of disk central density to cluster mean 
density, po,d{Rd) / pc- The heat input from the spheroid is 
negligible beyond 1 kpc for disk clusters and has therefore 
been ignored to provide scaling freedom. Surviving clusters 
are bounded by disruption at low density (small x)a.nA evap- 
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Table 4. Halo clusters: fraction of remaining mass after 10 Gyr as a function of orbital inclination 

10 5 M Q 



i = 


15° 


30° 


60° 


90° 


5 kpc 
10 kpc 


0.02 
0.49 


0.0 

0.59 


0.09 
0.71 


0.14 
0.73 




If 


8 M 






i = 
5 kpc 
10 kpc 


15° 
0.50 
0.58 


30° 
0.41 
0.69 


60° 
0.44 
0.80 


90° 
0.49 
0.80 


K = 1.0 and W = 


= 5 







limit set by the spheroid for Rd > 1 kpc but still strongly 
influenced by disk heating. Evolution is weakly dependent 
on the oscillation amplitude, with the maximum effect oc- 
curring between 2zo and 5zo- We exploit this property below 
to remove the vertical dimension from the phase space grid 
used to construct the distribution of evolved disk clusters. 

3.1.5 Concentration dependence 

The preceding discussion is based on clusters of a single 
concentration. To examine the concentration dependence, 
we consider the evolution of disk clusters of varying con- 
centration for a range of mass and at values of po.d(Rd) / pc 
which bracket the range of maximum cluster lifetime. Fig- 
ure ^| indicates that evaporation dominates at high concen- 
tration and high density while disruption dominates at low 
concentration and low density. Clusters show similar trends 
for increasing eccentricity, with evaporation becoming more 
important due to the increasing densities of tidally limited 
clusters at fixed Galactocentric radius. 




Figure 1. Remaining mass at 10 Gyr for disk clusters on circular 
orbits executing 5zo oscillations as a function of the ratio of disk 
central density to cluster central density (Wo = 5 profile). Solid 
contours show remaining masses in the range 3.0 < logM c < 
4.5 with AlogAf c = 0.5. Dotted contours show values of \ as 
labeled in the top right panel. Initial masses are given at the top 
of each panel. Lower density clusters tend to disrupt due to tidal 
heating while higher density clusters tend to evaporate due to 
rapid rates of relaxation. The density required for tidal limitation 
in the spheroid is too low to allow survival against disk heating. 

oration at high density (large %). Clusters with 10 5 Mq can 
only survive for Rd > 7 kpc. For higher mass, the evap- 
oration boundary moves to higher density because of the 
longer relaxation time scale. However, the disruption bound- 
ary remains roughly constant at density higher than the tidal 
limit. 

3.1.4 Dependence on oscillation height 

Table ^ shows the dependence of disk cluster evolution 
on oscillation height for clusters on circular orbits with 
log po.d(Rd)/ Pc = —0.37. This value is well within the tidal 



3.2 Internal properties 

To illustrate some basic trends in the evolution of internal 
cluster properties, we examine traits of the 10 6 Mq disk clus- 
ters with apogalactica at 8 kpc. Figure ^ compares evolution 
in projected profiles for cases ranging from tidal disruption 
to relaxation-dominated core contraction. The tidally dom- 
inated k — 1.0 and k — 0.9 clusters shows marked depar- 
ture from the initial profiles, developing a steeper, roughly 
power-law decline. The limiting radii remain near the ex- 
pected tidal boundary which moves inward due to the mass 
loss. For k = 0.7, heating is strong enough to produce de- 
viation from the initial King profile, although the central 
evolution is relatively unaffected. For « = 0.5, heating is so 
weak that the outer profile remains fixed while the central 
regions undergo gravothermal contraction. The tidally in- 
fluenced profiles also show mild concavity in the fall-off, a 
feature which resmbles the observed profiles given by Grill- 
mair et al. (1995,1996), who interpreted their observations 
as tidal tails. However, the feature evident here is not an 
unbound tidal tail but a bound halo region which has been 
partially cleared through orbital resonances. 

The evolution of the profile of the mass spectral index 
is shown in Figure W. Mass segregation occurs in every case, 
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Table 5. Disk clusters: fraction of remaining mass after 10 Gyr as a function of disk oscillation height 
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Figure 2. Remaining mass after 10 Gyr as a function of initial 
concentration and radius in the disk for clusters on circular or- 
bits with indicated masses and mean densities. Solid contours 
show remaining masses in the range 3.0 < logM c < 4.5 with 
AlogM c = 0.5. Evaporation/disruption isochrones at 5 Gyr (dot- 
ted) and 15 Gyr (dashed) are also shown. Clusters with 10 5 Mq 
evaporate for R^ < 6 kpc at high mean density (upper left) and 
evaporate or disrupt for R^ < 7 kpc at low mean density (lower 
left). Clusters with 10 6 Mq disrupt for R^ < 7 kpc and very low 
concentration at high mean density (upper right) and disrupt over 
a large range in radius and concentration for low mean density 
(lower right). 



Figure 3. The evolution of surface density profiles in four 10 6 Mq 
clusters on indicated orbits with apogalactica at 8 kpc. Solid lines 
show initial profiles; dotted lines show profiles after 5 Gyr; dashed 
lines show profiles after 10 Gyr except for re = 1.0 whose profile 
is shown at 7.5 Gyr, just prior to disruption. In the strong tidal 
cases (re = 1.0, re = 0.9), halos are truncated and profiles develop 
a steeper fall off than the initial profile. In the re = 0.7 case, no 
expansion occurs due to the near balance between heating and 
relaxation, but the outer profile evolves due to tidal heating. The 
re = 0.5 cluster undergoes negligible tidal heating, evidenced by 
the static halo profile, while relaxation leads to increasing central 
densities. 



regardless of the strength of tidal heating. Most low-mass 
stars evaporate while the k = 1.0 cluster disrupts as in- 
dicated by the strong flattening of /3(R) at all radii. The 
other cases show flattening of the spectrum in the core and 
steepening in the halo with differences that increase with ec- 
centricity. The increasing differences result from the shorter 
evolutionary time scales at high eccentricity for orbits with 
equal apocenter. Aside from differences in time scale, the 
evolution of the mass spectral index does not depend signif- 
icantly on orbit. The spectral index remains approximately 
constant in time near the initial half- mass radius of the clus- 
ter. However, the half-mass radius is relatively constant only 
in the most eccentric cases and undergoes considerable evo- 
lution whre tidal effects are strong. 



3.3 Characteristics of the present-day cluster 
population 

3. 3. 1 Distribution of cluster masses 

The mass spectrum of observed clusters is not well-described 
by a single power-law index over the mass range considered 
here. Harris & Pudritz (1994) find a change in slope near 
10 5 Mq . This trend is evident in Table [] which compares 
fits to the mass spectrum using three different models: a 
single power law, a two-component power law and a Gaus- 
sian magnitude distribution. The single power law shows a 
relatively flat spectrum in agreement with Harris & Pudritz 
(1994). The two-component power law shows a fairly steep 
dependence for M > 2.2 x 10 Mq and a nearly flat spectrum 
for masses below that. The Gaussian magnitude distribution 
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Figure 4. The line-of-sight mass spectral index as a function 
of projected cluster radius. Solid lines show the initial value, 
(3 = 2.35; dotted lines show the dependence after 5 Gyr; dashed 
lines show the dependence after 10 Gyr, except for k = 1.0 which 
shows (3{R) after 7.5 Gyr, just prior to disruption. Uncertainties 
are plotted in only one case to show the typical size of 1 — a error 
bars. The index (3{R) shows the effect of mass segregation in each 
case. The increased difference between core and halo indices with 
eccentricity results from the more rapid evolutionary timescale at 
high eccentricity for orbits with equal apocenter. 

peaks at 2.7 x 10 5 M Q , consistent with the two-component 
power law. 

Likelihood ratio tests show that the single power law 
can be rejected in favor of both the Gaussian magnitude 
distribution and two-component model at better than 99% 
confidence. The Gaussian and two-component power law can 
be discriminated with only 40% confidence. However, be- 
cause the data is so sparse, we adopt the single power law 
as the simplest model which provides a tenable description 
of the overall cluster mass distribution. 



3.3.2 Orbital dynamics and spatial distribution 

Figure ^ shows the inferred model parameters for the Mestel 
sphere both with and without the tangential velocities. The 
comparison shows the improved constraints that the addi- 
tional information provides. With the radial velocity data 
set, the tangential anisotropy of the Mestel sphere fit is un- 
constrained; in the augmented data set, the confidence con- 
tours close within q = 1.0 and r\ = 5.0 at 99% confidence. 
The contours are centered about isotropy and rule out strong 
anisotropy. We will use the augmented data set throughout 
the remainder of this paper. 

The estimated value of the anisotropy radius of the 
Eddington sphere, R a = 20kpc (Table hj), also indicates 
an isotropic distribution. A likelihood ratio test weakly fa- 
vors the Eddington sphere over the Mestel sphere with 75% 
confidence. The better fit results primarily because the Ed- 



Figure 5. Comparison of Mestel sphere fits to data with (solid) 
and without (dashed) available three-space velocities. The addi- 
tional information in the augmented data set provides some con- 
straint on the degree of tangential anisotropy compared to the 
pure radial velocity data. The uncertainty still dominates, but we 
can rule out strongly anisotropic distributions. 



dington sphere has a central core and a steep decline for 
R > 20kpc. However, as shown in Paper II, cluster distribu- 
tions probably follow steep power-law profiles initially and 
develop cores at later times through dynamical evolution. 
Both models can provide nearly coreless power-law density 
profiles, but the Eddington sphere becomes extremely radi- 
ally biased in this case. The Mestel sphere, by contrast, can 
have arbitrary orbital anisotropy for a given spatial profile. 
Therefore, considering the relative isotropy of the present- 
day population, the Mestel sphere provides a more realistic 
description of an initial cluster population. To investigate 
initial conditions below, we use only the Mestel sphere to 
model the spherical portion of the cluster distribution. 

The two-component analysis implies that 96% of the 
clusters are in a spherical component, according to the es- 
timated value of F (Table |^) . The inferred sphericity of the 
distribution is in agreement with previous analyses of the full 
cluster system (Frenk & White 1980; Thomas 1989; Chernoff 
& Djorgovski 1989). However, a likelihood ratio test rejects 
both the Mestel and Eddington spheres in favor of the two- 
component model at better than 99% confidence, indicating 
that the purely spherical models inadequately represent the 
dynamics of the full system. The preference stems from the 
strong tangential anisotropy attributed to the disk compo- 
nent. 

Strong tangential anisotropy is expected in a disk com- 
ponent due to the presence of the rapidly rotating system of 
high-metallicity disk clusters (Armandroff 1989). However, 
our estimate implies that only 4 clusters in the sample are as- 
sociated with the disk, while recent determinations indicate 
the presence of nearly 25 metal-rich disk clusters (Arman- 
droff 1993). This apparent contradiction may result from the 
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Table 6. Models of present-day mass spectrum 
Model «i 



' ai 



a 2 



T a2 M cut (M e ) 



power-law 
two-component 
power-law 



0.80 0.05 

0.034 0.15 1.62 0.15 



2.2 x 10 5 



Model 



a V 



(TV 



log L 



gaussian mag. 



-7.56 0.11 1.30 0.10 



-251.03 



Table 7. Comparison of fits to Mestel and Eddington spheres. 

Model rj a v q 



log L 



-260.56 
2.6 x 10 4 -248.38 



Mestel 2.89 0.25 

Eddington 2.54 0.10 



0.04 0.11 



20.1 



4.1 



-3544.4 
-3542.3 



obscuration of low-latitude clusters a nd th e greater evolu- 
tionary rate of clusters near the disk (§ 



3.1.2) 



both of which 

lead to a deficit of halo clusters near the disk. As a result, 
the 'hole' in the spherical component is filled by 'borrowing' 
the isotropic portion of the disk system which is kinemati- 
cally well-matched to the halo system. This leaves a residual 
system of disk clusters with very high angular momentum. 



3.4 Initial conditions of the cluster population 

Evolution of the pure Mestel sphere produces a profile which 
is shallower at present than in the past (Figure ]aj due to the 
more rapid evolution at small galactocentric radii. The ini- 
tial velocity distribution is tightly constrained to isotropic. 

Table H compares the best-fit parameters of the spher- 
ical model with the two-component model. The spheroid in 
the two-component model has a slightly shallower decline 
but is considerably more radially biased because higher an- 
gular momentum orbits are mainly associated with the disk. 
The uncertainties in the disk parameters are reduced be- 
cause the present-day population fraction (1-F) is larger 
than in the present-day model. The current population is 
estimated to number about 16 clusters. This is still some- 
what below the expected size, but again may be the result 
of imposing spherical symmetry on the evolved halo distri- 
bution. A likelihood ratio test rejects the purely spherical 
model with 97.5% confidence. 

The analysis predicts that the velocity distributions of 
both halo and disk components become more tangentially 
biased with time due to the more rapid evaporation of clus- 
ters on eccentric orbits. The initial orbit distribution of the 
halo component has fairly strong radial bias, with approxi- 
mately 41% of its kinetic energy in radial motion, while the 
disk component has a nearly isotropic initial orbit distribu- 
tion. In addition, the initial halo cluster density distribution 
has power law index of r\ — 2q = 3.38, while the disk cluster 
density distribution has power law index rjd — 2qd = 2.25. 

Figure |^ compares the cumulative distribution of clus- 
ters in our data sample with the evolved profile of the two- 
component model along with the separate contributions of 
the disk and sphere. The model matches the data fairly well. 
At small radii our models overestimate the expected number 



0.5 



-0.5 




Figure 6. Comparison of 50%,90%,95%,99% confidence lev- 
els in estimated initial parameters (solid) with levels in esti- 
mated present-day parameters superimposed from previous figure 
(dashed). Constraints are stronger on the initial distribution and 
indicate initial isotropy. The slope of the initial density distri- 
bution — (tj + 2q) is approximately r~ 3,35 while the present-day 
slope is best described as r~ 2 - 95 . 



of clusters. However, the KS confidence that the observed 
and model distribution differ is only 58%. 

Using the initial conditions from the two-component 
model, we derive the estimated initial distribution of clus- 
ters (Figure ^ . The total initial population has roughly 250 
clusters, with 200 in the spheroid and 50 in the disk. Ap- 
proximately 43% of the initial population remains. 
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Table 8. Present-day two-component model 

V 1 a Id Vz a d F log L 

estimate 2.75 -0.035 0.76 17.0 7.8 0.37 1.04 0.96 -3495.00 

uncertainty 0.24 0.12 0.05 7.08 3.34 0.16 0.13 tg'°t 



Table 9. Comparison of initial conditions in 2-component and spherical models 

spherical 

V <? a Vd <ld a d F log L 

estimate 3.35 -0.018 0.95 - -3486.94 
uncertainty 0.32 0.16 0.06 - 



2-component 

V <1 a 7)d qd a d F log L 

estimate 2.82 -0.28 0.96 2.38 0.066 1.15 0.89 -3478.79 

uncertainty 0.27 0.03 0.07 0.53 0.10 0.23 0.08 




R (kpc) R (kpc) 



Figure 7. Comparison of the evolved two-component model with 
the observed cumulative distribution of clusters with R < 30 kpc 
in the present sample. The total number of clusters is 108. The 
spherical component is estimated to have 92 clusters and the disk 
16. 



Figure 8. Comparison of initial (solid) and evolved (dashed- from 
previous figure) two-component models. The total initial number 
of clusters is estimated to be about 250 with 200 in the spheroid 
and 50 in the disk. 



4 DISCUSSION 

Before significant mass loss through stellar evolution, young 
clusters typically have strongly concentrated profiles (e.g. 
Elson, Fall & Freeman 1987). During this phase, concentra- 
tion is reduced as mass loss dominates gravothermal contrac- 
tion, driving the expansion of individual clusters (Chernoff & 
Weinberg 1990). When stellar evolution mass loss subsides — 
and our model calculations begin — the resulting population 
will consist of clust ers w ith a range of profiles. The cal- 
culations shown in §3.1.5 indicate, however, that evapora- 
tion times do not strongly depend on initial concentration. 
High concentration clusters undergo weaker tidal heating 



but maintain correspondingly larger rates of relaxation, so 
that evaporation times, for given mass, are approximately 
independent of concentration. The weak d ependence of evap- 
oration time on density is similar (§ 3.1.3[ ). The results of our 
calculations, therefore, do not depend significantly on the 
distribution of initial concentration. 

Estimates of initial population size are derived assum- 
ing that an initially spherical cluster distribution remains 
spherical. However, clusters on orbits confined near the disk 
evolve more rapidly than those on high inclination orbits, 
leaving an axisymmetric distribution with minimum density 
about the midplane of the disk. Estimates of initial halo 
cluster population size should not change significantly when 
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accounting for additional cluster loss at low inclination be- 
cause only ~ 10% more depletion of halo clusters is required 
to account for the remaining disk clusters. However, the ob- 
served disk cluster population is about 60% larger than our 
estimate of 16. Assuming that our estimate for the initial 
population size scales accordingly, we estimate that the disk 
population may have had roughly 80 members initially. 

The choice of the mass spectrum of clusters also in- 
fluences estimates of initial population size. While use of 
a single power law is necessitated by sample size, the two- 
component power law clearly provides the best fit. To gauge 
the importance of this choice, we use the two-component fit 
and the results of Paper II as a guide. We estimate that the 
spectral index may decrease by 0.2 for M > M cut and by 
0.3 for M < M cu t since the high mass range is very similar 
to that considered in Paper II while more rapid evolution 
occurs in the low mass range. This implies initial spectral 
indices of ai — 0.3 and 0:2 = 1.8 and 49% of the population 
with M < M C ut- In the single index fit with a = 0.95, 47% of 
the population has M < M cu t- Since most of these clusters 
evaporate within a Hubble time and since they dominate the 
depleted population, estimates of initial size will not change 
substantially when using a two-index model. 

Both the spatial and kinematic distributions of halo 
globular clusters differ from those of halo field stars. Harris 
(1976) and Zinn (1985) found that the spatial distribution 
of the full cluster system goes as 7?~ 3 ' 5 for 4 < R < 20 kpc, 
in good agreement with the profile of halo field stars (Free- 
man 1996). At smaller radii the profile becomes shallower 
and flattens into a core (Ostriker, Binney & Saha 1989). 
Similarly, while the distribution of halo field star orbits 
appears to have significant radial bias (Beers & Sommer- 
Larsen 1995), the distri butio n of halo cluster orbits appears 
to be nearly isotropic (§3.3.2). Our results from the present- 



day fits also exhibit the discrepancy between the spatial dis- 
tributions of halo clusters and halo field stars. The power- 
law profile derived using the Mestel sphere appears consid- 
erably flatter than R~ 3 S because we have included clusters 
with R < 4 kpc in the data set. The Eddington sphere best 
represents the distribution over the full radial range. It has a 
core, density p oc R~ 2 5 at small radius, a steep drop beyond 
20 kpc with density p oc R~ 4 ' 5 and the expected power-law 
decline with density p oc i? -3,5 at intermediate radii. 

Based on the estimates of initial conditions, we con- 
clude that cluster evolution can account for the difference 
between the spatial distributions of halo clusters and field 
stars. The estimated initial density profile p oc ^- 3 - 38 ± - 3 
in both spherical and two-component models over the full 
radial range of the data. This is in good agreement with the 
halo field star profile p oc ^- 3 - 29 ± - 24 recently determined by 
Sommer-Larsen & Zhen (1990) as well as the conventionally 
accepted value of R~ 3 ' 5 . 

Evolution also accounts for at least some differences 
in halo cluster and field star kinematics. As discussed in 



3.4, the predicted initial halo cluster population in the 2- 



component model has a strong radial bias which diminishes 
over time due to the preferential evaporation of clusters on 
eccentric orbits. Approximately 41% of the kinetic energy of 
this component is initially in radial motions. Compared to 
33% implied by the isotropy of the observed halo cluster dis- 
tribution, this is closer to the value of 54% derived for halo 



field stars by Beers & Sommer-Larsen (1995) and is roughly 
consistent with the value of 44% derived by Norris (1986). 

The predicted initial disk population, by contrast, is un- 
like the observed stellar thick disk. The estimates indicate 
that the present-day high angular momentum population de- 
veloped from a nearly isotropic initial distribution through 
the selective evaporation of clusters on eccentric orbits. Al- 
though the model underestimates the size of the disk popu- 
lation, any resulting kinematic bias predicts higher angular 
momentum in the system because lower angular momentum 
members cannot be distinguished from the halo clusters. Ev- 
idence for a flattened component with nearly random kine- 
matics may exist in the sample of halo clusters studied by 
Zinn (1993) and in the halo field star sample studied by 
Sommer-Larsen & Zhen (1990). 

Overall, the results of this paper point to a scenario 
which correlates the origin of halo field stars and clusters 
during Galaxy formation. The predicted disk cluster system 
is less clearly correlated with observed disk stellar popu- 
lations but may be related to an intermediate phase of the 
dissipative collapse which is thought to have given rise to the 
disk (Larson 1990; Zinn 1993). The results do not contra- 
dict merger-induced heating of an initially cold disk (Quinn, 
Hernquist & Fullager 1993), provided that initially circular 
disk cluster orbits can be isotropized by the accretion event. 

Finally, we emphasize that our predictions describe the 
initial population of clusters which evolve quasi-statically 
through relaxation and tidal heating. Since clusters prob- 
ably formed with a range of initial densities, some would 
have been prey to rapid tidal disruption (Paper I). How- 
ever, the significance of such initial conditions cannot be 
divined from the observed population because, excepting a 
few clusters such as Pal 5 (Cudworth 1993), the majority 
of observed clusters must have formed with initial densities 
which allowed quasi-static evolution and long-term survival. 
Consequently, the importance of tidal disruption in shaping 
the cluster population at an early epoch can be described 
only through models of the initial conditions of the cluster 
system in a proto-Galactic or cosmological context (Paper 

I). 



5 CONCLUSIONS 

The main conclusions of this work are as follows: 

(i) The disk dominates tidal heating on low-eccentricity 
orbits. The effect of the spheroid is negligible except for 
R < 1 kpc. 

(ii) Disk oscillations heat more efficiently than disk shock- 
ing due to better frequency matching in typical cluster po- 
tentials. 

(iii) Cluster evolution depends weakly on initial concen- 
tration and density due to the combined effects of tidal heat- 
ing and evaporation 

(iv) The evolution of disk clusters depends weakly on os- 
cillation height. 

(v) The loss of clusters on low-inclination orbits implies 
that the distribution of halo clusters is less dense near the 
disk. 

(vi) The evaporation of clusters on high-eccentricity or- 
bits leads to increasing tangential bias in the distribution of 
cluster orbits. 
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(vii) The estimated initial spatial distribution of halo 
clusters matches the present-day distribution of halo field 
stars. The estimated initial kinematic distribution is nearer 
to the observed radial anisotropy in the kinematics of halo 
field stars, having approximately 40% of its energy in radial 
motions. 

(viii) The estimated initial distribution of disk clusters 
does not match the kinematic distribution of the stellar disk. 
However, it is similar to the flattened halo samples studied 
by Zinn (1993) and Sommer-Larsen & Zhen (1990). 

(ix) The outer profiles of evolving clusters deviate from 
the initial profiles in cases of moderate to strong tidal heat- 
ing, leading to considerable evolution in the half-mass ra- 
dius. Profiles may also show density inflections which result 
from resonant clearing and which resemble the observed pro- 
files of Grillmair et al. (1995,1996). 

(x) Profiles of the mass spectral index become natter in 
the core and steeper in the halo due to mass segregation. 
The profiles evolve similarly on all orbits and differ at fixed 
times only through orbitally determined differences in evo- 
lutionary timescale. 
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APPENDIX A: HEATING RATE FOR DISK 
OSCILLATIONS 

Clusters confined to the disk oscillate about the midplane 
and are heated through resonance with the periodic tidal 
compression in the direction perpendicular to the disk plane. 
As discussed in Weinberg (1994), the tidal potential of the 
disk is well-approximated by the leading term of an expan- 
sion of the disk potential about the center of mass of the 
cluster: 
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Hi = 2irGp[Z(t)]; 



(Al) 



where we have substituted the density for the second deriva- 
tive of the disk potential using Poisson's equation. The quan- 
tity Z(t) is the cluster position as a function of time, and z 
refers to the position of a star relative to the center of the 
cluster. 

Expanding the z 2 factor in the action-angle Fourier se- 
ries as defined in Tremaine & Weinberg (1984), we obtain 



(A2) 



Substituting the resulting expression for H\ into equation 
(5) of Paper I, we derive the heating rate 

oo 

((E)) = -8n 4 p4L ]T s h0 (± + ±8i 2 o)(i-n) 2 \xY 2 \ 2 

1 — — ~x; 

oo 

J2 \a n \ 2 S(l-tl-nLu).(A3) 

n = — oo 

Choosing the vertical profile of the disk and the amplitude 
of vertical oscillations completely specifies the rate of energy 
input for given cluster profile. In the present work, we adopt 
a Gaussian vertical profile for the disk. Comparison with an 
exponential vertical profile reveals little difference in overall 
heating rate at any oscillation amplitude. 

APPENDIX B: COORDINATE SYSTEMS 

Our analysis uses three coordinate systems: heliocentric 
spherical coordinates, Galactocentric spherical coordinates 
and Galactocentric cylindrical coordinates. The first system 
is the natural observational reference frame while the second 
and third are the natural reference frames for the sphere and 
disk models, respectively. In the heliocentric frame, we use 
(r, I, b) or distance, Galactic longitude and Galactic latitude. 
In the Galactocentric spherical frame, we use the (R, "!>, O) 
to denote Galactocentric distance, colatitude and azimuth 
respectively. In the Galactocentric cylindrical frame, we use 
the (R d ,<&,Z) to denote Galactocentric radius in the disk, 
azimuth and height above the plane, respectively. 

Velocity components in a particular direction are de- 
noted with the corresponding subscript. In the heliocentric 
frame, (iv, Vi, vi) are the radial, azimuthal and latitudinal 
components, respectively. In the Galactocentric spherical 
frame, (vr, v&,vq) are the radial, azimuthal and latitudinal 
components, respectively. In the Galactocentric cylindrical 
frame, (vR d ,v<s>,vz) are the polar radial, azimuthal and ver- 
tical velocities, respectively. We also write the components 
in the equivalent inner product form so that, for example, 
the Heliocentric velocity components are (y ■ f", v ■ I, v ■ b). 

APPENDIX C: GENERALIZED MESTEL DISK 

The phase space distribution function for disk clusters in the 
thin disk approximation has the form 

dN 



dE d 8J 2 dE z 



f d (E d ,J 2 )g(E z ) 



(CI) 



where f d (E d , J z ) governs the distribution in the plane of 
the disk, g(E z ) governs the distribution perpendicular to 
the disk, E d denotes orbital energy in the disk, J z refers to 
the angular momentum along the Z-axis and E z denotes the 
energy of vertical osillations. 

The Mestel disk distribution is defined as 



fd{E d ,J 2 z ) = Ae- E « /rT *J 2qd , 



(C2) 



where a d is the isothermal velocity dispersion in the disk. 
The radial and azimuthal velocity dispersions are (j\ d = a 2 d 
and a% = (2q d + l)a d which implies q d — a%/2a 2 ld — i, so 
that — ^ < q d < oo . The quantity q d defines the anisotropy of 
the orbit distribution in the disk. For q d — — i, the distribu- 
tion has only radial orbits while as q d — > oo the distribution 
has only circular orbits (Henon 1973; Barnes et al. 1986). 
The isothermal vertical distribution is defined as 



g(E z ) = Be 



(C3) 



where a z is the isothermal vertical velocity dispersion. In 
the Milky Way, a 2 varies with radius in the disk because the 
scale height is constant while the midplane density varies. 
Assuming it to be fixed introduces some bias into the ex- 
pected vertical velocities at a given radius. However, this 
has no effect on our conclusions because we find that heat- 
ing by disk oscillations is nearly independent of oscillation 
amp litude (or, equivalently, velocity at the midplane; c.f. 

Integrating over va d , v<s>, vz and Z yields the surface 
density in the logarithmic potential 

= cV^(2* d ) q+1 r(q + \)R d ^- 2 ^ (C4) 

where we absorb all vertical integration constants into the 
factor C and define the parameter r] d 
2, the density goes as ln_Rd- 



2 ,Ja 2 d . Forr/ d -2q d = 



APPENDIX D: 
SPHERE 



GENERALIZED MESTEL 



The phase space distribution function for halo clusters has 
the form 
dN 



f(E, r 



(Dl) 



8E8J 2 

where E is the total energy and J is the total angular mo- 
mentum of a cluster. 

The Mestel sphere has the same form as the Mestel disk, 



I(e, r 



Ae- E/a J 2 \ 



(D2) 



but is a three-dimensional distribution so that the radial 
velocity dispersion a R — a and tangential velocity disper- 
sion crfn = 2(q + 1)<7 2 . Consequently, q = a 2 ^/2aj l — 1 and 
— 1 < q < oo. The quantity q defines the anisotropy of the 
orbit distribution. For q = — 1, the distribution has only of 
radial orbits while as q — » oo the distribution has only circu- 
lar orbits (Henon 1973; Barnes et al. 1986). Integrating over 
vr, ii$ and we gives the volume density in the logarithmic 
potential 



(D3) 



4^ = A-K 3/2 {2o 2 Y +s/2 Y{q + l)R-^~ 2 ^ 
a^R 

where we define r) = v 2 /a 2 . For 77 — 2q = 3, the density goes 
as lni?. 
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APPENDIX E: GENERALIZED EDDINGTON 
SPHERE 



The Eddington model has the distribution function 



f{E, J 2 



Ae 



-E/<r 2 



-J 2 /2R 2 cr 2 



(El) 



which implies that a 2 R — a 2 and cr\ — 2<r 2 /(l + R 2 /R 2 ). 
Integrating over velocities gives the volume density in the 
logarithmic potential 



d3R - A2n(2na ) 1 + R2/R z 
where we define r\ = v 2 /a 2 . 



(E2) 



(Aguilar, Hut & Ostriker 1988; Thomas 1989). For exam- 
ple, if we take a spherical model for the cluster distribution 
function f(E, J 2 ), the marginal probability of any given ob- 
servation is 



f(v r ,R)= I dv l( iv b f(-(vt +vt + vt) + 9{R), 

R 2 ((v-$) 2 + (v-e) 2 )), 



(F3) 



where 9(R) is the potential, and the tangential velocity com- 
ponents are written in inner product notation. 



APPENDIX F: LIKELIHOOD WITH 
INCOMPLETE DATA SETS 

An incomplete data set is defined as one in which some ob- 
servations are missing. For each cluster, the observations 
have varying completeness with respect to the full set of ob- 
servations required in a given model. There are several ap- 
proaches to deriving parameter estimates in this situation. 
In the present work, we adopt a likelihood-based estimation 
scheme^. 

If we denote the full phase-space vector for our mod- 
els (f,v,m) as Y and write Y = {Y obs ,Y mis ), where Y obs 
signifies the observed data and Y m i S signifies the missing 
data, then f(Y\0) = f(Y oba ,Y m i S \9) denotes the underlying 
probability of observing all quantities Y ts and Y m i S , where 
f(-\9) is governed by the parameters 9. Integrating the dis- 
tribution over each case of missing data gives the marginal 
probability density of Y b 3 : 



f(Y obs \9) = / f{y obs ,Y mi3 \9)dY m 



(Fl) 



For independent observations, we denote the marginal prob- 
ability for the i th cluster as fi = f(Yi tt >bs\Q) and write the 
likelihood function in the usual way as the joint probability 
of the observations given the model: 



(F2) 



Using L(9) to derive inferences concerning 9 requires 
that there are no selection effects leading to the system- 
atic absence of data for a particular class of observations. 
In practice, of course, we know that this is not the case 
for observations of globular clusters, where, for example, 
latitude-dependent extinction results in the absence of radial 
velocities. However, in the present analysis we make no at- 
tempt to derive any type of selection function. One approach 
which does not ignore selection effects is the Expectation- 
Maximization algorithm, an iterative procedure which pro- 
vides estimates for the model parameters as well as the miss- 
ing data (Little & Rubin 1987, ch. 7). 

The standard data set used in analyses of the spatial 
distribution and kinematics of the cluster system consists 
of cluster positions, masses and heliocentric radial velocities 



t The following discussion is based on the presentation of Little 
& Rubin (1987), chapter 5. 
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